// CVisImage.hpp
// Developed by Yang Yu (yuy@janelia.hhmi.org), Sept 27, 2012
// 

#ifndef __CVISIMG_HPP
#define __CVISIMG_HPP

#include "CVisImage.h"

//------------------------------------
/// nd image (XYZCT...) structure
//------------------------------------

template<class Tdata, class Tidx>
Y_IMAGE<Tdata, Tidx> :: Y_IMAGE()
{
    dk=IMAGE;
    dt=(DataType)sizeof(Tdata);
    dimt=UNKNOWNDIMTYPE;
    pImg=NULL;
    dims=NULL;
    centroid=NULL;
    totalplxs=0;
}

template<class Tdata, class Tidx>
Y_IMAGE<Tdata, Tidx> :: ~Y_IMAGE()
{
}

// clean
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: clean()
{
    y_del<Tdata>(pImg);
    y_del<Tidx>(dims);

    y_del<REAL>(centroid);

    return;
}

// set an image
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: setImage(Tdata *p, Tidx *sz, Tidx imgdim) //
{
    if(!p || !sz || imgdim<=0)
    {
        cout<<"Invalid inputs"<<endl;
        return;
    }

    pImg=p; dims=sz; dimt = (DimType)imgdim;

    totalplxs=1;
    for(Tidx i=0; i<dimt; i++) totalplxs*=dims[i];

    if(dimt==4) set4doffsets(); // 4D image

    refc = 0;

    return;
}

// init an image
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: initImage(Tidx *sz, Tidx imgdim) // init a zero image
{
    if(!sz || imgdim<=0)
    {
        cout<<"Invalid inputs"<<endl;
        return;
    }
    dimt = (DimType)imgdim;
    y_new<Tidx, Tidx>(dims, imgdim);

    totalplxs=1;
    for(Tidx i=0; i<imgdim; i++)
    {
        dims[i] = sz[i];
        totalplxs *= dims[i];
    }

    try
    {
        pImg=new Tdata [totalplxs];
        memset(pImg, 0, sizeof(Tdata)*totalplxs);
    }
    catch(...)
    {
        printf("Fail to allocate memory.\n");
        y_del<Tdata>(pImg);
        return;
    }

    if(dimt==4) set4doffsets(); // 4D image

    refc = 0;

    return;
}

// init a field
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: initField(Tidx *sz, Tidx imgdim) // init a zero field
{
    dk = FIELD;
    initImage(sz, imgdim);
    return;
}

// set 4d image offset
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: set4doffsets()
{
    if(dims)
    {
        of3 = dims[2]*dims[1]*dims[0];
        of2 = dims[1]*dims[0];
        of1 = dims[0];
    }
    else
    {
        cout<<"Fail to set offsets"<<endl;
    }

    return;
}

// get mass center
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: getMassCenter()
{
    if(pImg && dims)
    {
        y_new<REAL, Tidx>(centroid, (Tidx)dimt);

        REAL sum = 0.0;
        switch(dimt)
        {
        case D1D:

            // X

            break;
        case D2D:

            // XY

            break;
        case D3D:

            // XYZ

        case D4D:

            // XYZC C=1

            centroid[0] = 0.0;
            centroid[1] = 0.0;
            centroid[2] = 0.0;

            for(Tidx z=0; z<dims[2]; z++)
            {
                for(Tidx y=0; y<dims[1]; y++)
                {
                    for(Tidx x=0; x<dims[0]; x++)
                    {
                        REAL val = (REAL)(val4d(refc, z, y, x));

                        centroid[0] += val*(REAL)x;
                        centroid[1] += val*(REAL)y;
                        centroid[2] += val*(REAL)z;

                        sum += val;
                    }
                }
            }

            centroid[0] /= sum;
            centroid[1] /= sum;
            centroid[2] /= sum;

            break;
        case D5D:

            break;
        default:
            printf("Unsupported dimensions.\n");
        }
    }
}

// set reference channel
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: setRefChn(Tidx refcinput)
{
    refc = refcinput;
}

// zeros
template<class Tdata, class Tidx>
void Y_IMAGE<Tdata, Tidx> :: zeros()
{
    for(Tidx i=0; i<totalplxs; i++) pImg[i] = 0;
}

// voxel value of 4d image
template<class Tdata, class Tidx>
Tdata Y_IMAGE<Tdata, Tidx> :: val4d(Tidx z, Tidx y, Tidx x, Tidx v)
{
    return (pImg[ z*of3 + y*of2 + x*of1 + v ]);
}

// 1d index of 4d image
template<class Tdata, class Tidx>
Tidx Y_IMAGE<Tdata, Tidx> :: idx4d(Tidx z, Tidx y, Tidx x, Tidx v)
{
    return (z*of3 + y*of2 + x*of1 + v);
}


#endif //__CVISIMG_HPP
